library(phyloseq)
library(vegan)
library(dplyr)
library(reshape2)
library(here)
library(ggplot2)
library(DESeq2)
library(ggsci)
library(Biostrings)
library(seqinr)
library(janitor)
library(cowplot)
library(forcats)
library(tidyr)
library(DECIPHER)
library(phangorn)
library(ggrepel)
npg <- pal_d3("category10")(10)
theme_char <- function(base_size = 11, base_family = ""){
theme_bw() %+replace%
theme(axis.text = element_text(color = "Black"))
}
theme_set(theme_char())
set.seed(12)
ps <- readRDS(here("marcos_flea", "ps_trial.rds"))
sample_data <- read.csv(here("marcos_flea", "Sample_Info.csv"))
rownames(sample_data) <- sample_data$ID
sample_data(ps) <- sample_data
psprop <- ps %>%
transform_sample_counts(function(otu) otu/sum(otu))
multi <- genefilter_sample(ps, filterfun_sample(function(x) x > 0), A = 2)
multips <- prune_taxa(multi, ps)
multips.glom <- tax_glom(psprop, taxrank = rank_names(ps)[6])
tax_table <- as.data.frame(tax_table(ps))
tax_table$OTU <- rownames(tax_table)
seq <- refseq(ps)
ps
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 483 taxa and 10 samples ]
## sample_data() Sample Data: [ 10 samples by 22 sample variables ]
## tax_table() Taxonomy Table: [ 483 taxa by 7 taxonomic ranks ]
## refseq() DNAStringSet: [ 483 reference sequences ]
ps.glom <- tax_glom(ps, taxrank = rank_names(ps)[6])
ps.melt <- psmelt(ps)
psp.glom <- tax_glom(psprop, taxrank = rank_names(ps)[6])
pspg.melt <- psmelt(psp.glom)
psp.melt <- psprop %>% psmelt
otu_table <- as.data.frame(otu_table(ps))
otu_sums <- as.data.frame(colSums(otu_table))
colnames(otu_sums) <- "Abundance"
otu_sums$OTU <- rownames(otu_sums)
tax_table <- left_join(tax_table, otu_sums, by = "OTU")
wolb.melt <- psp.melt %>% filter(Genus == "Wolbachia")
wolbotu <- unique(wolb.melt$OTU)
wolb.melt$Big3 <- ifelse(wolb.melt$OTU == "ASV1", "ASV1", "Other")
wolb.melt$Big3 <- ifelse(wolb.melt$OTU == "ASV2", "ASV2", wolb.melt$Big3)
wolb.melt$Big3 <- ifelse(wolb.melt$OTU == "ASV3", "ASV3", wolb.melt$Big3)
wolb_otu_tab <- otu_table %>% select(wolbotu)
ps.wolb <- prune_taxa(wolbotu, ps)
bart.melt <- ps.melt %>% filter(Genus == "Bartonella")
bartotu <- unique(bart.melt$OTU)
allfam <- names(sort(taxa_sums(psp.glom), decreasing=TRUE))
top20 <- names(sort(taxa_sums(psp.glom), decreasing=TRUE))[1:20]
ps.top20 <- transform_sample_counts(psp.glom, function(OTU) OTU/sum(OTU))
ps.top20 <- prune_taxa(top20, ps.top20)
plot_bar(ps.top20, x="ID", fill="Family")+ facet_wrap(~Time, scales="free_x")
ggsave(here("marcos_flea", "Image", "top20.png"), plot = last_plot(), width = 9, height = 6)
tax_table %>% dplyr::filter(Genus == "Wolbachia") %>% dplyr::mutate(sum = sum(Abundance))
## Kingdom Phylum Class Order Family
## 1 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 2 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 3 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 4 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 5 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 6 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 7 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 8 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 9 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 10 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 11 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 12 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 13 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 14 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 15 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 16 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 17 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 18 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 19 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 20 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 21 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 22 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 23 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 24 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 25 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 26 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 27 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 28 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 29 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 30 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 31 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 32 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 33 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 34 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## Genus Species OTU Abundance sum
## 1 Wolbachia <NA> ASV1 1336358 3159186
## 2 Wolbachia <NA> ASV2 1331247 3159186
## 3 Wolbachia <NA> ASV3 392785 3159186
## 4 Wolbachia <NA> ASV9 30562 3159186
## 5 Wolbachia <NA> ASV10 30379 3159186
## 6 Wolbachia <NA> ASV22 17357 3159186
## 7 Wolbachia <NA> ASV41 9373 3159186
## 8 Wolbachia <NA> ASV131 2705 3159186
## 9 Wolbachia <NA> ASV134 2673 3159186
## 10 Wolbachia <NA> ASV232 1324 3159186
## 11 Wolbachia <NA> ASV233 1319 3159186
## 12 Wolbachia <NA> ASV390 413 3159186
## 13 Wolbachia <NA> ASV450 298 3159186
## 14 Wolbachia <NA> ASV468 280 3159186
## 15 Wolbachia <NA> ASV502 246 3159186
## 16 Wolbachia <NA> ASV538 212 3159186
## 17 Wolbachia <NA> ASV564 191 3159186
## 18 Wolbachia <NA> ASV596 167 3159186
## 19 Wolbachia <NA> ASV599 164 3159186
## 20 Wolbachia <NA> ASV609 159 3159186
## 21 Wolbachia <NA> ASV617 157 3159186
## 22 Wolbachia <NA> ASV679 128 3159186
## 23 Wolbachia <NA> ASV714 113 3159186
## 24 Wolbachia <NA> ASV789 91 3159186
## 25 Wolbachia <NA> ASV844 71 3159186
## 26 Wolbachia <NA> ASV850 69 3159186
## 27 Wolbachia <NA> ASV864 62 3159186
## 28 Wolbachia <NA> ASV951 46 3159186
## 29 Wolbachia <NA> ASV960 45 3159186
## 30 Wolbachia <NA> ASV962 45 3159186
## 31 Wolbachia <NA> ASV1011 38 3159186
## 32 Wolbachia <NA> ASV1013 38 3159186
## 33 Wolbachia <NA> ASV1021 37 3159186
## 34 Wolbachia <NA> ASV1039 34 3159186
3159186/sum(otu_table)
## [1] 0.9666248
fam_replace <- tax_table
fam_replace$Family <- replace_na(fam_replace$Family, "Unassigned")
fam_freq <- fam_replace %>%
# mutate(Family = fct_infreq(Family)) %>%
ggplot(aes(x = Family))+
geom_bar()+
coord_flip()+
labs(y = "Number of ASVs")
fam_abund <- fam_replace %>% group_by(Family) %>% mutate(fam_freq = sum(Abundance)) %>%
distinct(Family, .keep_all = TRUE) %>%
ggplot(aes(x = Family, y = fam_freq))+
geom_col()+
coord_flip()+
labs(y = "Abundance", x = "")+theme(axis.text.y = element_blank())
plot_grid(fam_freq, fam_abund, rel_widths = c(1, 0.5))
ggsave(here("marcos_flea", "Image", "Fig2.eps"), width = 8, height = 12.5)
gen_glom <- fam_replace %>% group_by(Genus) %>% mutate(fam_freq = sum(Abundance)) %>%
distinct(Family, .keep_all = TRUE)
fam_glom <- fam_replace %>% group_by(Family) %>% mutate(fam_freq = sum(Abundance)) %>% mutate(fam_asvs = n()) %>% distinct(Family, .keep_all = TRUE)
tax_glom <- as.data.frame(tax_table(ps.glom))
tax_glom$OTU <- rownames(tax_glom)
prop_table <- as.data.frame(otu_table(ps.glom))
prop_sums <- as.data.frame(colSums(prop_table))
colnames(prop_sums) <- "Abundance"
prop_sums$OTU <- rownames(prop_sums)
otu_prop <- prop_sums %>% mutate(Proportion = (round(Abundance/sum(Abundance)*100, 4)))
taxa_prop <- left_join(tax_glom, otu_prop, by = "OTU")
print_genus <- data.frame(Family = taxa_prop$Family,
Genus = taxa_prop$Genus,
Proportion = taxa_prop$Proportion)
write.csv(print_genus, here("marcos_flea", "genus_prop.csv"))
invsimpson <- as.data.frame(vegan::diversity(otu_table, index = "invsimpson"))
colnames(invsimpson) <- "invsimpson"
invsimpson$ID <- rownames(invsimpson)
shannon <- as.data.frame(vegan::diversity(otu_table, index = "shannon"))
colnames(shannon) <- "shannon"
shannon$ID <- rownames(shannon)
spprich <- as.data.frame(specnumber(otu_table)) #calc species richness
colnames(spprich) <- "spprich"
spprich$ID <- rownames(spprich)
even <- as.data.frame(shannon[,1]/log(spprich[,1])) #calc Pielou's evenness
colnames(even) <- "even"
even$ID <- rownames(spprich)
diversity <- left_join(invsimpson, shannon, by = "ID")
diversity <- left_join(diversity, spprich, by = "ID")
diversity <- left_join(diversity, even, by = "ID")
diversity <- left_join(diversity, sample_data, by = "ID")
diversity$Cat_Bart <- ifelse(diversity$Cat_Bart == "Bartonella Naive", "Bartonella Uninfected", diversity$Cat_Bart)
diversity$group <- ifelse(diversity$ID == "50UF", "50UF", diversity$Cat)
diversity$group <- ifelse(diversity$ID == "47UF", "47UF", diversity$group)
diversity$Cat_Bart <- ifelse(diversity$ID == "47UF", "Unfed Flea", diversity$Cat_Bart)
diversity$Cat_Bart <- ifelse(diversity$ID == "50UF", "Unfed Flea", diversity$Cat_Bart)
shannon.gg <- diversity %>%
ggplot(aes(x = as.factor(Time), y = shannon, color = Cat_Bart, group = group))+
geom_point(alpha = 0.8)+geom_line()+
labs(x = "Time (days)", y = "Shannon Index")+
theme(axis.text.x = element_text(), legend.position = "none")+
scale_color_manual(values = c("red", "black", "grey50"))
invsimpgg <- diversity %>%
ggplot(aes(x = as.factor(Time), y = invsimpson, color = Cat_Bart, group = group))+
geom_point(alpha = 0.8)+geom_line()+
labs(x = "Time (days)", y = "Inverse Simpson Index")+
theme(axis.text.x = element_text(), legend.position = "none")+
scale_color_manual(values = c("red", "black", "grey50"))
spprichgg <- diversity %>%
ggplot(aes(x = as.factor(Time), y = spprich, color = Cat_Bart, group = group))+
geom_point(alpha = 0.8)+geom_line()+
labs(x = "Time (days)", y = "Species Richness")+
theme(axis.text.x = element_text(), legend.position = "none")+
scale_color_manual(values = c("red", "black", "grey50"))
evennessgg <- diversity %>%
ggplot(aes(x = as.factor(Time), y =even, color = Cat_Bart, group = group))+
geom_point(alpha = 0.8)+geom_line()+
labs(x = "Time (days)", y = "Pielou's Evenness")+
theme(axis.text.x = element_text(), legend.position = "none")+
scale_color_manual(values = c("red", "black", "grey50"))
bart.type <- c(
expression(paste(italic("Bartonella"), " Infected")),
expression(paste(italic("Bartonella"), " Uninfected")),
"Unfed Flea"
)
fill.tit <- expression(paste("Cat ", italic("Bartonella"), " Status"))
legend_plot <- diversity %>%
ggplot(aes(x = as.factor(Time), y =even, color = Cat_Bart, group = group))+
geom_point()+geom_line()+
labs(color = fill.tit)+
scale_color_manual(values = c("red", "black", "grey50"), labels = bart.type)+
theme(legend.text.align = 0)
legend <- get_legend(legend_plot)
ugh <- plot_grid(spprichgg, evennessgg, shannon.gg, legend, nrow = 1, labels = c("A", "B", "C"))
save_plot(plot = ugh, here("marcos_flea", "Image", "Fig3.jpeg"), base_width = 8, base_height = 3.5)
ugh
bhshan <- diversity %>% filter(Cat_Bart == "Bartonella Infected")
noshan <- diversity %>% filter(Cat_Bart == "Bartonella Uninfected")
t.test(bhshan$shannon, noshan$shannon, paired = FALSE)
##
## Welch Two Sample t-test
##
## data: bhshan$shannon and noshan$shannon
## t = 1.403, df = 3.6663, p-value = 0.2394
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.3082353 0.8941888
## sample estimates:
## mean of x mean of y
## 1.543527 1.250551
bhshan <- diversity %>% filter(Cat_Bart == "Bartonella Infected") %>% filter(Time == "1")
noshan <- diversity %>% filter(Cat_Bart == "Bartonella Uninfected")%>% filter(Time == "1")
t.test(bhshan$shannon, noshan$shannon)
##
## Welch Two Sample t-test
##
## data: bhshan$shannon and noshan$shannon
## t = 3.0834, df = 1.9442, p-value = 0.0943
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2370271 1.3280859
## sample estimates:
## mean of x mean of y
## 1.872497 1.326967
d1shan <- diversity %>% filter(Time == "1")
d9shan <- diversity %>% filter(Time == "9")
t.test(d1shan$shannon, d9shan$shannon, paired = T)
##
## Paired t-test
##
## data: d1shan$shannon and d9shan$shannon
## t = 2.5038, df = 3, p-value = 0.08742
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.1098851 0.9206561
## sample estimates:
## mean difference
## 0.4053855
bray <- vegdist(otu_table, "bray")
bray <- melt(as.matrix(bray), varnames = c("row", "col"))
ggplot(bray, aes(x = row, y = col, fill = value))+
geom_tile()
NMDS = metaMDS(otu_table, k = 2)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.04473556
## Run 1 stress 0.09181492
## Run 2 stress 0.1015291
## Run 3 stress 0.0506377
## Run 4 stress 0.09411252
## Run 5 stress 0.05063741
## Run 6 stress 0.05063741
## Run 7 stress 0.04473556
## ... Procrustes: rmse 5.726807e-06 max resid 1.358413e-05
## ... Similar to previous best
## Run 8 stress 0.0506373
## Run 9 stress 0.04473556
## ... Procrustes: rmse 3.930524e-06 max resid 5.979458e-06
## ... Similar to previous best
## Run 10 stress 0.04473556
## ... Procrustes: rmse 2.150821e-05 max resid 4.993466e-05
## ... Similar to previous best
## Run 11 stress 0.04473556
## ... Procrustes: rmse 5.359859e-06 max resid 9.428021e-06
## ... Similar to previous best
## Run 12 stress 0.04473556
## ... Procrustes: rmse 2.865025e-06 max resid 6.84322e-06
## ... Similar to previous best
## Run 13 stress 0.04473556
## ... Procrustes: rmse 1.628114e-05 max resid 3.832866e-05
## ... Similar to previous best
## Run 14 stress 0.04473556
## ... Procrustes: rmse 2.123653e-06 max resid 3.700492e-06
## ... Similar to previous best
## Run 15 stress 0.04473556
## ... Procrustes: rmse 2.346207e-05 max resid 5.489858e-05
## ... Similar to previous best
## Run 16 stress 0.05063741
## Run 17 stress 0.08960182
## Run 18 stress 0.05063743
## Run 19 stress 0.09411241
## Run 20 stress 0.09181554
## *** Best solution repeated 8 times
sample.scores <- as.data.frame(scores(NMDS)[1])
sample.scores$ID <- rownames(sample.scores)
sample.scores <- left_join(sample.scores, sample_data, by = "ID")
asv.scores <- as.data.frame(scores(NMDS)[2])
sample.scores$xnudge <- "0.05"
sample.scores$ynudge <- "0.05"
sample.scores$xnudge <- ifelse(sample.scores$ID == "3711-24hr", 0.10, sample.scores$xnudge)
sample.scores$Cat_Bart <- ifelse(sample.scores$Cat_Bart == "Bartonella Uninfected", "Bartonella Uninfected", sample.scores$Cat_Bart)
sample.scores$ID <- ifelse(sample.scores$ID == "50UF", "UF50", sample.scores$ID)
sample.scores$ID <- ifelse(sample.scores$ID == "47UF", "UF47", sample.scores$ID)
sample.scores$Cat_Bart <- ifelse(sample.scores$ID == "UF47", "Unfed Flea", sample.scores$Cat_Bart)
sample.scores$Cat_Bart <- ifelse(sample.scores$ID == "UF50", "Unfed Flea", sample.scores$Cat_Bart)
bart.type <- c(
expression(paste(italic("Bartonella"), " Infected")),
expression(paste(italic("Bartonella"), " Uninfected")),
"Unfed Flea"
)
fill.tit <- expression(paste("Cat ", italic("Bartonella"), " Status"))
ggplot(sample.scores, aes(x = sites.NMDS1, y = sites.NMDS2, color = Cat_Bart))+
geom_point()+
theme_char()+
geom_text(aes(label = ID), nudge_x = 0.23)+
labs(color = fill.tit)+
guides(color = guide_legend(override.aes=list(shape = 15, size = 5)))+
scale_color_manual(values = c("red", "black", "grey45"), labels = bart.type)+
theme(legend.text.align = 0)
ggsave(here("marcos_flea", "Image", "Fig4.eps"), plot = last_plot(), width = 7, height = 4.5)
bart.melt %>%
ggplot(aes(x = ID, y = Abundance, fill = Cat_Bart))+
geom_col(color = "black")+theme_char()+
facet_grid(OTU~Cat_Bart, scales = "free_x")+
theme(axis.text.x = element_text(angle = 45,hjust = 1), legend.position = "none")+
labs(y = "Abundance")+
geom_text(aes(label = Abundance), vjust = -0.30)
ggsave(here("marcos_flea", "Image", "bart.png"), plot = last_plot(), width = 9, height = 6)
wolbseq <- seq[wolbotu]
writeXStringSet(wolbseq, here("marcos_flea", "sequences", "wolbseq.fasta"))
pspg.melt %>% filter(Genus == "Wolbachia") %>%
ggplot(aes(y = Abundance, x = ID, fill = as.factor(Cat_Bart)))+
geom_col(color = "black")+theme_char()+
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
labs(y = "Wolbachia Proportion", color = "Cat Bartonella")+
facet_grid(~Time, scales = "free_x")+
scale_fill_manual(values = c("Pink", "Black"))
wolb.melt %>% mutate(Cat = factor(Cat, levels = c("Unfed", "3363", "3508", "3320", "3711"))) %>%
ggplot(aes(x = ID, y = Abundance, fill = Big3))+
geom_col(color = "black")+theme_char()+
facet_wrap(~Cat, scales = "free_x", nrow = 1)+
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
labs(y = "Proportion of Total Reads")
ggsave(here("marcos_flea", "Image", "wolbbig3.png"), plot = last_plot(), width = 9, height = 6)
for(i in 1:10){
graph <- wolb.melt %>% filter(OTU == wolbotu[i]) %>% filter(Abundance > 0) %>%
ggplot(aes(x = as.factor(Time), y = Abundance, color = Cat, group = Cat))+
geom_point()+geom_line()+
facet_wrap(~Cat_Bart, scales = "free_x")+theme_char()+
labs(x = "Time (days)", y = "Proportion", title = paste(wolbotu[i], "Proportion"))
print(graph)
}
wolb.interest <- c("ASV1", "ASV2", "ASV3")
wolb.melt$group <- ifelse(wolb.melt$Sample == "50UF", "50UF", wolb.melt$Cat)
wolb.melt$group <- ifelse(wolb.melt$Sample == "47UF", "47UF", wolb.melt$group)
wolb.melt$group <- ifelse(wolb.melt$Sample == "47UF", "47UF", wolb.melt$group)
wolb.melt$Cat_Bart <- ifelse(wolb.melt$Sample == "47UF", "Unfed Flea", wolb.melt$Cat_Bart)
wolb.melt$Cat_Bart <- ifelse(wolb.melt$Sample == "50UF", "Unfed Flea", wolb.melt$Cat_Bart)
wolb.melt$Cat_Bart <- ifelse(wolb.melt$Cat_Bart == "Bartonella Naive", "Bartonella Uninfected", wolb.melt$Cat_Bart)
bart.type <- c(
expression(paste(italic("Bartonella"), " Infected")),
expression(paste(italic("Bartonella"), " Uninfected")),
"Unfed Flea"
)
fill.tit <- expression(paste("Cat ", italic("Bartonella"), " Status"))
wolb.melt %>% filter(OTU %in% wolb.interest) %>% mutate() %>%
mutate(OTU = factor(OTU, levels = c(wolb.interest))) %>%
mutate(Cat_Bart = factor(Cat_Bart, levels = c("Bartonella Infected", "Bartonella Uninfected", "Unfed Flea"))) %>%
ggplot(aes(x = as.factor(Time), y = Abundance, color = Cat_Bart, group = group))+
geom_line()+geom_point()+
facet_wrap(~OTU)+theme_char()+
labs(x = "Time (days)", y = "Proportion", color = fill.tit)+
scale_color_manual(values = c("red","black", "grey65"), labels = bart.type)+
theme(legend.position = "right", legend.text.align = 0)
ggsave(here("marcos_flea", "Image", "Fig5.tiff"), plot = last_plot(), width = 7, height = 3.5)
asv.interest <- c("ASV1", "ASV2", "ASV3", "ASV9", "ASV10", "ASV352")
psp.melt$group <- ifelse(psp.melt$Sample == "50UF", "50UF", psp.melt$Cat)
psp.melt$group <- ifelse(psp.melt$Sample == "47UF", "47UF", psp.melt$group)
psp.melt$group <- ifelse(psp.melt$Sample == "47UF", "47UF", psp.melt$group)
psp.melt$Cat_Bart <- ifelse(psp.melt$Sample == "47UF", "Unfed Flea", psp.melt$Cat_Bart)
psp.melt$Cat_Bart <- ifelse(psp.melt$Sample == "50UF", "Unfed Flea", psp.melt$Cat_Bart)
psp.melt %>% filter(OTU %in% asv.interest) %>% mutate() %>%
mutate(OTU = factor(OTU, levels = c(asv.interest))) %>%
mutate(Cat_Bart = factor(Cat_Bart, levels = c("Bartonella Naive", "Bartonella Infected", "Unfed Flea"))) %>%
ggplot(aes(x = as.factor(Time), y = Abundance, color = Cat_Bart, group = group))+
geom_line()+geom_point()+
facet_wrap(~OTU, scales = "free_y")+theme_char()+
labs(x = "Time (days)", y = "Proportion", color = "Host Bartonella Status")+
scale_color_manual(values = c("black", "red", "grey69"))+
theme(legend.position = "bottom")
hm <-prune_taxa(names(sort(taxa_sums(ps),TRUE)[1:300]), ps)
plot_heatmap(hm, sample.label="ID", na.value = "white")
Cat Bart w/o nonfed
multinounfed <- subset_samples(multips, Fed == "Yes")
deseq <- phyloseq_to_deseq2(multinounfed, ~ Cat_Bart)
## converting counts to integer mode
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
deseq <- DESeq(deseq, test = "Wald", fitType = "parametric")
## estimating size factors
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## Note: levels of factors in the design contain characters other than
## letters, numbers, '_' and '.'. It is recommended (but not required) to use
## only letters, numbers, and delimiters '_' or '.', as these are safe characters
## for column names in R. [This is a message, not a warning or an error]
## final dispersion estimates
## fitting model and testing
res <- results(deseq, cooksCutoff = FALSE)
alpha = 0.01
res = cbind(as(res, "data.frame"), as(tax_table(multips)[rownames(res), ], "matrix"))
wolb.dseq <- res %>% filter(Genus == "Wolbachia") %>% filter(padj < 0.05)
wolb.dseq
## [1] baseMean log2FoldChange lfcSE stat pvalue
## [6] padj Kingdom Phylum Class Order
## [11] Family Genus Species
## <0 rows> (or 0-length row.names)
sigtab = res[which(res$padj < alpha),]
head(sigtab)
## baseMean log2FoldChange lfcSE stat pvalue padj
## ASV288 127.53884 -24.79779 3.098097 -8.004200 1.202456e-15 1.394849e-13
## ASV336 77.19327 -24.11462 3.383905 -7.126269 1.031263e-12 3.987551e-11
## ASV377 59.77222 -23.76011 3.384060 -7.021185 2.199948e-12 5.103880e-11
## ASV391 76.36428 -23.88336 3.383911 -7.057917 1.690174e-12 4.901504e-11
## ASV402 31.38034 -22.89174 3.384681 -6.763339 1.348472e-11 2.214814e-10
## ASV432 24.36797 -21.79153 3.385056 -6.437570 1.214016e-10 1.280235e-09
## Kingdom Phylum Class Order
## ASV288 Bacteria Firmicutes Clostridia Lachnospirales
## ASV336 Bacteria Firmicutes Clostridia Lachnospirales
## ASV377 Bacteria Firmicutes Clostridia Lachnospirales
## ASV391 Bacteria Firmicutes Clostridia Lachnospirales
## ASV402 Bacteria Firmicutes Clostridia Lachnospirales
## ASV432 Bacteria Actinobacteriota Coriobacteriia Coriobacteriales
## Family Genus Species
## ASV288 Lachnospiraceae Oribacterium asaccharolyticum
## ASV336 Lachnospiraceae [Ruminococcus] torques group <NA>
## ASV377 Lachnospiraceae <NA> <NA>
## ASV391 Lachnospiraceae Blautia stercoris
## ASV402 Lachnospiraceae [Eubacterium] fissicatena group <NA>
## ASV432 Eggerthellaceae Eggerthella lenta
# Phylum order
x = tapply(sigtab$log2FoldChange, sigtab$Phylum, function(x) max(x))
x = sort(x, TRUE)
sigtab$Phylum = factor(as.character(sigtab$Phylum), levels=names(x))
# Genus order
x = tapply(sigtab$log2FoldChange, sigtab$Genus, function(x) max(x))
x = sort(x, TRUE)
sigtab$Genus = factor(as.character(sigtab$Genus), levels=names(x))
sigtab$OTU <- rownames(sigtab)
sigtab$Analysis <- c("Cat Bartonella Status")
ggplot(sigtab, aes(x=Genus, y=log2FoldChange, color=Phylum)) +
geom_point(size=2) +
theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))+coord_flip()+
labs(title = "Cat Bartonella Status")
Time w/0 unfed
deseq.time <- phyloseq_to_deseq2(multinounfed, ~ Time)
## converting counts to integer mode
## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() function
deseq.time <- DESeq(deseq.time, test = "Wald", fitType = "parametric")
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
res.time <- results(deseq.time, cooksCutoff = FALSE)
alpha = 0.05
res.time = cbind(as(res.time, "data.frame"), as(tax_table(multinounfed)[rownames(res.time), ], "matrix"))
wolb.time <- res.time %>% filter(Genus == "Wolbachia") %>% filter(padj < 0.05)
wolb.time
## [1] baseMean log2FoldChange lfcSE stat pvalue
## [6] padj Kingdom Phylum Class Order
## [11] Family Genus Species
## <0 rows> (or 0-length row.names)
sigtab.time = res.time[which(res.time$padj < alpha),]
sigtab.time$OTU <- rownames(sigtab.time)
head(sigtab.time)
## baseMean log2FoldChange lfcSE stat pvalue padj
## ASV193 145.90173 -3.563588 0.3815665 -9.339363 9.691577e-21 1.124223e-18
## ASV279 68.48759 -3.428953 0.3789772 -9.047915 1.457247e-19 8.452035e-18
## ASV314 75.77528 -3.422823 0.3828236 -8.940993 3.856871e-19 1.223514e-17
## ASV336 77.19327 2.468604 0.4229935 5.836034 5.345796e-09 3.875702e-08
## ASV367 68.10692 -3.424755 0.3834652 -8.931072 4.219012e-19 1.223514e-17
## ASV377 59.77222 2.985442 0.4230134 7.057559 1.694531e-12 1.310438e-11
## Kingdom Phylum Class Order Family
## ASV193 Bacteria Firmicutes Clostridia Oscillospirales Ruminococcaceae
## ASV279 Bacteria Firmicutes Clostridia Lachnospirales Lachnospiraceae
## ASV314 Bacteria Firmicutes Clostridia Lachnospirales Lachnospiraceae
## ASV336 Bacteria Firmicutes Clostridia Lachnospirales Lachnospiraceae
## ASV367 Bacteria Firmicutes Clostridia Lachnospirales Lachnospiraceae
## ASV377 Bacteria Firmicutes Clostridia Lachnospirales Lachnospiraceae
## Genus Species OTU
## ASV193 Incertae Sedis <NA> ASV193
## ASV279 Lachnoclostridium <NA> ASV279
## ASV314 <NA> <NA> ASV314
## ASV336 [Ruminococcus] torques group <NA> ASV336
## ASV367 Lachnospiraceae NK4A136 group <NA> ASV367
## ASV377 <NA> <NA> ASV377
# Phylum order
x.time = tapply(sigtab.time$log2FoldChange, sigtab.time$Phylum, function(x) max(x))
x.time = sort(x.time, TRUE)
sigtab.time$Phylum = factor(as.character(sigtab.time$Phylum), levels=names(x.time))
# Genus order
x.time = tapply(sigtab.time$log2FoldChange, sigtab.time$Genus, function(x) max(x.time))
x.time = sort(x.time, TRUE)
sigtab.time$Genus = factor(as.character(sigtab.time$Genus), levels=names(x.time))
sigtab.time$Analysis <- c("Time")
sigtab.time$OTU <- rownames(sigtab.time)
ggplot(sigtab.time, aes(x = 1, y = OTU, fill = log2FoldChange))+
geom_tile()+
theme(axis.text.y = element_text(size = 6), axis.text.x = element_blank())+
facet_wrap(~Phylum, scales = "free_y")
ggplot(sigtab.time, aes(x=Genus, y=log2FoldChange, color=Phylum)) +
geom_point( alpha = 0.8) +
theme(axis.text.x = element_text(),
legend.position = "bottom")+coord_flip()
Fed Status
deseq.fed <- phyloseq_to_deseq2(multips, ~ Fed)
## converting counts to integer mode
deseq.fed <- DESeq(deseq.fed, test = "Wald", fitType = "parametric")
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 43 genes
## -- DESeq argument 'minReplicatesForReplace' = 7
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
res.fed <- results(deseq.fed, cooksCutoff = FALSE)
res.fed = cbind(as(res.fed, "data.frame"), as(tax_table(multips)[rownames(res.fed), ], "matrix"))
wolb.fed <- res.fed %>% filter(Genus == "Wolbachia") %>% filter(padj < 0.05)
wolb.fed
## baseMean log2FoldChange lfcSE stat pvalue padj
## ASV390 12.35088 18.66602 3.88657 4.802698 1.565417e-06 3.506535e-05
## Kingdom Phylum Class Order
## ASV390 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales
## Family Genus Species
## ASV390 Anaplasmataceae Wolbachia <NA>
sigtab.fed = res.fed[which(res.fed$padj < alpha),]
alpha = 0.01
sigtab.fed$OTU <- rownames(sigtab.fed)
head(sigtab.fed)
## baseMean log2FoldChange lfcSE stat pvalue padj
## ASV285 23.91761 19.57942 3.885561 5.039019 4.679249e-07 1.903797e-05
## ASV314 24.58057 19.61182 3.885533 5.047394 4.478771e-07 1.903797e-05
## ASV367 14.15733 18.85508 3.886304 4.851674 1.224237e-06 3.427863e-05
## ASV390 12.35088 18.66602 3.886570 4.802698 1.565417e-06 3.506535e-05
## ASV400 14.36153 17.76983 3.886279 4.572453 4.820462e-06 7.712740e-05
## ASV403 23.06325 19.51557 3.885603 5.022533 5.099456e-07 1.903797e-05
## Kingdom Phylum Class Order
## ASV285 Bacteria Firmicutes Clostridia Oscillospirales
## ASV314 Bacteria Firmicutes Clostridia Lachnospirales
## ASV367 Bacteria Firmicutes Clostridia Lachnospirales
## ASV390 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales
## ASV400 Bacteria Firmicutes Clostridia Lachnospirales
## ASV403 Bacteria Acidobacteriota Acidobacteriae Acidobacteriales
## Family Genus Species OTU
## ASV285 Ruminococcaceae Incertae Sedis <NA> ASV285
## ASV314 Lachnospiraceae <NA> <NA> ASV314
## ASV367 Lachnospiraceae Lachnospiraceae NK4A136 group <NA> ASV367
## ASV390 Anaplasmataceae Wolbachia <NA> ASV390
## ASV400 Lachnospiraceae GCA-900066575 <NA> ASV400
## ASV403 <NA> <NA> <NA> ASV403
# Phylum order
x.fed = tapply(sigtab.fed$log2FoldChange, sigtab.fed$Phylum, function(x) max(x))
x.fed = sort(x.fed, TRUE)
sigtab.fed$Phylum = factor(as.character(sigtab.fed$Phylum), levels=names(x.fed))
# Genus order
x.fed = tapply(sigtab.fed$log2FoldChange, sigtab.fed$Genus, function(x) max(x.fed))
x.fed = sort(x.fed, TRUE)
sigtab.fed$Genus = factor(as.character(sigtab.fed$Genus), levels=names(x.fed))
sigtab.fed$Analysis <- c("Fed Status")
ggplot(sigtab.fed, aes(x = 1, y = OTU, fill = log2FoldChange))+
geom_tile()+
theme(axis.text.y = element_text(size = 6), axis.text.x = element_blank())+
facet_wrap(~Phylum, scales = "free_y")
ggplot(sigtab.fed, aes(x=Genus, y=log2FoldChange, color=Phylum)) +
geom_point( alpha = 0.8) +
theme(axis.text.x = element_text(),
legend.position = "right")+coord_flip()+
labs(title = "Fed Status")
All Together
alldseq <- bind_rows(sigtab, sigtab.fed, sigtab.time)
alldseq <- alldseq[complete.cases(alldseq$Genus),]
ggplot(alldseq, aes(x=OTU, y=Analysis, fill = log2FoldChange)) +
geom_tile()+
theme(axis.text.x = element_text(),
legend.position = "right",
axis.text.y = element_text(size = 10))+coord_flip()
ggplot(alldseq, aes(x=Genus, y=log2FoldChange, fill = Phylum)) +
geom_col(position = "dodge")+
theme(axis.text.x = element_text(),
legend.position = "bottom",
axis.text.y = element_text(size = 10))+coord_flip()+
facet_wrap(~Analysis)
alldseq %>% mutate(Analysis = factor(Analysis, levels = c("Fed Status", "Cat Bartonella Status", "Time"))) %>%
mutate(OTU = paste0(Genus, " (",OTU, ")")) %>%
ggplot(aes(x=OTU, y=log2FoldChange)) +
geom_segment(aes(x = OTU, y = 0, xend = OTU, yend = log2FoldChange))+
geom_point(size = 2, aes(color = Order))+
theme(axis.text.x = element_text(),
legend.position = "bottom",
axis.text.y = element_text(size = 10))+
coord_flip()+
facet_wrap(~Analysis)
ggsave(here("marcos_flea", "Image", "alldeseq.png"), plot = last_plot(), width = 8, height = 5.5)
psp.melt %>% filter(OTU %in% alldseq$OTU) %>% filter(Abundance > 0) %>%
mutate(OTU = paste0(Genus, " (",OTU, ")")) %>%
ggplot(aes(x = as.factor(Time), y = Abundance, color = Cat_Bart))+
geom_point()+
facet_wrap(~OTU, nrow = 5)+
theme(axis.text.x = element_text(angle = 90, hjust = 1), legend.position = "bottom", strip.text = element_text(size = 11))+
labs(x = "Time (Days)", y = "Proportion of Reads")+
scale_color_manual(values = c("Red", "Black", "grey"))
ggsave(here("marcos_flea", "Image", "individualdeseq.png"), plot = last_plot(), width = 13, height = 7)
psprop %>% psmelt %>% filter(Genus == "Blautia") %>% filter(Abundance > 0) %>%
ggplot(aes(y = Abundance, x = as.factor(Time), fill = Cat_Bart))+
geom_col(aes(color = Cat_Bart))+geom_point(color = "Black", alpha = 0.6)+
theme(axis.text.x = element_text(angle = 90, hjust = 1), legend.position = "right")+
labs(x = "Time (Days)", y = "Proportion of Reads", title = "Blautia Abundance")